Homework 11

第十一次作業
b94611001 歐陽太閒
本週(5/24)有來上課

一,凸輪資料分析:

本次作業的凸輪之資料如下
升程:
加速度運動:0-100°滯留後開始提升(rise)到5cm,並且在200-260°間滯留,而後進入返程.
返程:
加速度運動:260-320°間為返程至0cm,滯留至360°
刻度:
10°

以上之狀況可做如下分析(加速度運動)θ為輸入角度,β=區間

升程上半段位移:y(θ)=2h(θ/β)²
升程下半段位移:y(θ)=h[1-2(1-θ/β)²]
返程上半段位移:y(θ)=h[1-2(1-θ/β)²]
返程下半段位移:y(θ)=2h(1-θ/β)²

因此我們可以利用老師講義中的parabol_cam.m來進行此凸輪在不同角度時的位置,速度,加速度對角度的關係:

function [y, yy, yyy]=parabol_cam(phi, phi_in,beta_range,direct,travel,rpm)
各項參數分別的代表意義為:
phi=當前輸入角度
phi_in=升降程起始角度
beta_range=本程所經之角度
direct=方向,1為升程,-1為返程
travel=衝程
rpm=轉速

因此可以對於凸輪的運動進行分析,在此選擇返程也為一加速度運動,當然也可以選擇其他運動.

程式架構:
1.程式資訊
2.建立向量
3.計算升程引用parabol_cam
4.計算停滯
5.計算返程引用parabol_cam
6.計算停滯
7.向量合成
8.繪圖

故撰寫程式如下:

function [ss, vv, aa]=acccam(rstart,rend,rstart2,rend2,c,stroke)
%作者 b94611001 歐陽太閒
%rstart=升程起始角度,rend=升程結束角度,rstart2=降程開始角度,rend2=降程結束角度,c=分割度數,stroke=衝程
%for example:acccam(100,200,260,320,10,5)

title('作業11 作者:歐陽太閒 b94611001'); %資訊
xlabel('x-axis degree');
ylabel('y-axis cm');
axis equal;
itmr=rstart/c; %進行角度處理
itme=rend/c;
itmr2=rstart2/c;
itme2=rend2/c;
hold on; %多層繪圖

theta1=linspace(0,0,itme); %建立一0向量

for q=1:itmr %先填入升程前的角度
theta1(q)=q*c
end
for j=itmr:length(theta1); %計算升程
theta1(j)=j*c; %填入角度
[ss1(j), vv1(j), aa1(j)]=parabol_cam(theta1(j),rstart,rend-rstart,1,stroke,0); %引用parabol_cam
end

theta2=linspace(rend,rstart2,itmr2-itme+1); %計算滯留
ss2=linspace(ss1(itme),ss1(itme),itmr2-itme+1); %延續升程資料
vv2=linspace(vv1(itme),vv1(itme),itmr2-itme+1);
aa2=linspace(0,0,itmr2-itme+1);

theta3=rstart2:c:rend2 %計算返程
for k=1:(length(theta3))
[ss3(k), vv3(k), aa3(k)]=parabol_cam(theta3(k),rstart2,rend2-rstart2,-1,5,0) %引用parabol_cam
end;

theta4=linspace(rend2,360,(360/c)-itme2+1); %計算滯留
ss4=linspace(ss3(length(theta3)),ss3(length(theta3)),(360/c)-itme2+1);
vv4=linspace(vv3(length(theta3)),vv3(length(theta3)),(360/c)-itme2+1);
aa4=linspace(0,0,(360/c)-itme2+1);

theta=[theta1 theta2 theta3 theta4];
ss=[ss1 ss2 ss3 ss4]; %合併向量
vv=[vv1 vv2 vv3 vv4];
aa=[aa1 aa2 aa3 aa4];
h1=plot(theta,ss,'b-',theta,vv,'k-',theta,aa,'r-') %繪製
legend('位移','速度','加速度',3) %標註


因此我們可以輸入以下指令

>>acccam(100,200,260,320,10,5)

結果如圖


只顯示位移時:
結果如圖



當然我們也可以使用老師的plot_dwell.m來進行繪製,但plotdwell只適用於末段為回程者,但我們也可以測試之.
升程:
加速度運動:0-100°滯留後開始提升(rise)到5cm,並且在200-260°間滯留,而後進入返程.
返程:
加速度運動:260-360°間為返程至0cm
刻度:
10°

可以輸入以下指令

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

結果如圖

二,凸輪工作曲線

在此作業中,需繪製第一題中的凸輪之工作曲線半徑為15公分,順時鐘方向旋轉,梢形垂直接觸從動件,長10cm
因此本題可以老師的

function [x,y]=pincam(cth,r0,s,e,L,range,pattern,cw)

為基礎,由於選擇的範圍dwell.m無法處理,故無法直接使用pincam.m,因此改寫一程式pincam_acc,其中引用了上一題中的acccam.m計算結果,而

可得此凸輪之工作曲線.

程式架構:
1.程式資訊
2.處理偏置量
2.引入acccam計算資料(位移)
3.計算轉軸並存入向量
4.繪製延伸線及梢
5.繪圖

程式碼如下:

function [x,y]=pincam_acc(cth,r0,s,e,L,rstart,rend,rstart2,rend2,cw)
%作者 b94611001 歐陽太閒
%rstart=升程起始角度,rend=升程結束角度,rstart2=降程開始角度,rend2=降程結束角度,cth=範圍,s=衝程,r0=凸輪基圓半徑
%e=偏置量,L:從動件(銷型)長度,cw:凸輪轉動方向(+1,-1)
%for example:pincam_acc([0:10:360],15,5,0,10,100,200,260,320,+1)
title('作業11 作者:歐陽太閒 b94611001'); %資訊
xlabel('x-axis cm');
ylabel('y-axis cm');
th=cth*pi/180;
s0=sqrt(r0*r0-e*e); %含偏置量時計算
[ss,vv,aa]=acccam(rstart,rend,rstart2,rend2,10,s) %引入acccam計算資料
for i=1:length(cth) %引入範圍內的位移,並以軸旋轉
t=th(i)*cw;
A=[cos(t) -sin(t);sin(t) cos(t)];%以軸轉
ym=ss(i);
x0=s0+ym;
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','blue') %繪製梢
end
hold on;
plot([0 x],[0 y],'ro',x,y,'k-') %繪圖
rb = linspace(0, 2*pi, 301);
plot(r0*cos(rb), r0*sin(rb),'r-'); %畫基圓
axis equal;


輸入pincam_acc([0:10:360],15,5,0,10,100,200,260,320,-1)

>>pincam_acc([0:10:360],15,5,0,10,100,200,260,320,-1)

結果如圖


顯示基圓時:



而假設有偏置量e=5時

>>pincam_acc([0:10:360],15,5,5,10,100,200,260,320,-1)

結果如圖


本題若選擇為
升程:
加速度運動:0-100°滯留後開始提升(rise)到5cm,並且在200-260°間滯留,而後進入返程.
返程:
加速度運動:260-360°間為返程至0cm

符合dwell.m的設計,則可使用pincam.m

輸入
pincam([0:10:360],15,5,0,10,[100 200 260],[2 2],-1)
結果如圖


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

與剛剛撰寫的pincam_acc([0:10:360],15,5,0,10,100,200,260,360,-1)結果相同

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

結果如圖

三,凸輪之動畫

在本程式中,嘗試利用rotate(把手,[1 0 0],角度)這個命令,可以大幅減少程式碼,並且嘗試更簡潔的動畫結構.

程式結構如下:
1.程式資訊
2.處理偏置量
2.引入acccam計算資料(位移)
3.計算轉軸並存入向量
4.進行動畫繪製
a.旋轉及平移
b.暫停


因此撰寫以下程式:

function [x,y]=pincammov(cth,r0,s,e,L,rstart,rend,rstart2,rend2,cw)
%作者 b94611001 歐陽太閒
%rstart=升程起始角度,rend=升程結束角度,rstart2=降程開始角度,rend2=降程結束角度,cth=範圍,s=衝程,r0=凸輪基圓半徑
%e=偏置量,L:從動件(銷型)長度,cw:凸輪轉動方向(+1,-1)
%for example:pincam_acc([0:10:360],15,5,0,10,100,200,260,320,+1)
title('作業11 作者:歐陽太閒 b94611001'); %資訊
xlabel('x-axis cm');
ylabel('y-axis cm');
th=cth*pi/180;
s0=sqrt(r0*r0-e*e); %含偏置量時計算
[ss,vv,aa]=acccam(rstart,rend,rstart2,rend2,10,s) %引入acccam計算資料

for i=1:length(cth) %引入範圍內的位移,並以軸旋轉
t=th(i)*cw;
A=[cos(t) -sin(t);sin(t) cos(t)];%以軸轉
ym=ss(i);
x0=s0+ym;
Sx=[0 x0 x0+L;e e e];
X=A\Sx;
x(i)=X(1,2);y(i)=X(2,2);
end

h1=plot([0 x],[0 y],'ro',x,y,'k-') %標點
axis equal;
axis([-25 25 -25 25]);

for j=1:100 %進行旋轉
q=37-rem(j,37);
ys=r0+ss(q);
%clf;
h2=line([0 0],[ys ys+L],'linewidth',3,'color','blue'); %繪製梢
rotate(h1,[0 0 1],(-10));
pause(0.1);
h2=line([0 0],[ys ys+L],'linewidth',3,'color','white');
pause(0.05)
end

輸入pincam_acc([0:10:360],15,5,0,10,100,200,260,320,-1)

>>pincam_acc([0:10:360],15,5,0,10,100,200,260,320,-1)

結果如圖


影片


附錄:acccam(100,200,260,320,10,5) 運算結果

ss =
Columns 1 through 7
0 0 0 0 0 0 0
Columns 8 through 14
0 0 0 0.1000 0.4000 0.9000 1.6000
Columns 15 through 21
2.5000 3.4000 4.1000 4.6000 4.9000 5.0000 5.0000
Columns 22 through 28
5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000
Columns 29 through 35
4.7222 3.8889 2.5000 1.1111 0.2778 0.0000 0.0000
Columns 36 through 39
0.0000 0.0000 0.0000 0.0000
vv =
Columns 1 through 7
0 0 0 0 0 0 0
Columns 8 through 14
0 0 11.4592 12.6051 13.7510 14.8969 16.0428
Columns 15 through 21
5.7296 4.5837 3.4377 2.2918 1.1459 0 0
Columns 22 through 28
0 0 0 0 0 0 -82.7606
Columns 29 through 35
-85.9437 -89.1268 -9.5493 -6.3662 -3.1831 -0.0000 -0.0000
Columns 36 through 39
-0.0000 -0.0000 -0.0000 -0.0000
aa =
Columns 1 through 7
0 0 0 0 0 0 0
Columns 8 through 14
0 0 6.5656 6.5656 6.5656 6.5656 6.5656
Columns 15 through 21
-6.5656 -6.5656 -6.5656 -6.5656 -6.5656 -6.5656 0
Columns 22 through 28
0 0 0 0 0 0 -18.2378
Columns 29 through 35
-18.2378 -18.2378 18.2378 18.2378 18.2378 18.2378 0
Columns 36 through 39
0 0 0 0

Homework 10

作業10
本週(5/17)有來上課
b94611001 歐陽太閒

一,速度與加速度

本次作業需要分析速度與加速度,故可參考網路講義第三章和第六章,講義中有三種作法化成向量:
1.複數法
R=r*e^(jθ)
2.尤拉法
R=r*(cosθ+jsinθ)
以下採用較易在直角座標系分析的尤拉法進行分析

1.1位置

當一桿R以某特定點M等角速度ω迴轉時,其端點P之速度方向
我們可考慮一桿R,有水平角度θ,
θ=ω*t + 0.5*α*t^2 +a
其尤拉向量式如下

R = r*cosθ + j*r*sinθ

因此我們可以寫出下列程式碼:

M=[0 0];
x=M(1)+ r*cos(theta) %求出R位置的x,y分量
y=M(2)+ r*sin(theta)
plot([M(1) x],[M(2) y]) %畫出

例如:角度=60,r=5

d2r=pi/180;
theta=60;
r=5;
M=[0 0];
x=M(1)+ r*cos(theta); %求出R位置的x,y分量
y=M(2)+ r*sin(theta);
plot([M(1) x],[M(2) y]); %畫出位置
%或者
%linkshape([M(1) M(2)],[x y]);

或者使用1.4的函數

>>mec_alys(5,60,0,0,[0 0])

x =
2.5000

y =
4.3301


結果如圖

1.2速度(以下微分假設r長不變)

令α=0,a=0而求桿端速度則需將向量微分(ds/dt),即可得
v = -r*ω*sinθ + j*r*ω*cosθ
微積分厲害之處在此顯現,速度的水平分量其中即已包含逆時針旋轉時和桿水平向量相反的符號,不必進行額外的處理,
只需要平移即可,因此其方向為:

u=-sinθ + jcosθ

因此我們可以寫出以下的程式碼:
omg為角速度,
例如:角度=60,r=5,角速度=1

theta=60*d2r
r=5
omg=1
vx= -r*sin(theta)*omg %將位置向量微分算出速度向量
vy= r*cos(theta)*omg
u= -sin(theta)+i*cos(theta) %速度方向
plot([M(1) x],[M(2) y])
%或者
%linkshape([M(1) M(2)],[x y]);
plot([x x+vx],[y y+vy],'color','r') %平移後繪出
text(x+vx,y+vy,['u=' num2str(-1*sin(theta),'%6.1f') '+i' num2str(cos(theta),'%6.1f')]) %顯示


或者使用1.4的函數

>>mec_alys(5,60,1,0,[0 0])
x =
2.5000
y =
4.3301
vx =
-4.3301
vy =
2.5000
u =
-0.8660 + 0.5000i %方向

1.3加速度

我們要求加速度則需將速度微分(dv/dt)並整理後可得

a=(-rωcosθ+rαsinθ)+j(-rωsinθ+rαcosθ)

因此可寫出以下程式

ax=-1*r*(omg^2)*cos(theta) + r*alpha*sin(theta)
ay=-1*r*(omg^2)*sin(theta) + r*alpha*cos(theta)
plot([x+vx x+vx+ax],[y+vy y+vy+ay],'color','r') %平移後繪出
if ax<0.001; %防止除零
ax=0.001;
end;
dir=atan(ay/ax); %求方向
au=cos(dir) + i*sin(dir);

例如:角度=60,r=5,角速度=1,角加速度=1.2

>>mec_alys(5,60,1,1.2,[0 0])
x =
2.5000
y =
4.3301
vx =
-4.3301
vy =
2.5000
u =
-0.8660 + 0.5000i %角速度方向
ax =
2.6962
ay =
-1.3301
au =
0.8968 - 0.4424i %角加速度方向


結果如圖

1.4化成函數

因此我們得到

R = r*cosθ + j*r*sinθ
v = -r*ω*sinθ + j*r*ω*cosθ
a = (-rω2cosθ+rαsinθ)+j(-rω2sinθ+rαcosθ)

而其方向為

uv=-sinθ + jcosθ
av= cos(atan((-rω2sinθ+rαcosθ)/-rω2cosθ+rαsinθ))+sin(atan((-rω2sinθ+rαcosθ)/-rω2cosθ+rαsinθ))

故我們可以合併以上分析(1.1,1.2,1.3)寫出以下函數:


程式架構:
1.資訊
2.求位置
求速度,計算方向
求加速度,計算方向
3.繪製

function mec_alys(l,theta,omega,alpha,Ori)
%作者 b94611001 歐陽太閒
%l=桿長,theta=角度,omega=角速度,ori=原點
title('作業10 作者:歐陽太閒 b94611001'); %資訊
xlabel('x-axis cm');
ylabel('y-axis cm');
axis equal;
d2g=pi/180; %角徑轉換
r=l; %桿長
M=Ori; %原點
omg=omega; %角速度
hold on;
theta=theta*d2g;
costh=cos(theta);
sinth=sin(theta);
x=M(1)+ r*costh
y=M(2)+ r*sinth
%求速度
nsinth=-1*sinth;
vx= r*nsinth*omg
vy= r*costh*omg
u=nsinth+i*costh %速度方向
%加速度
ax=-1*r*(omg^2)*costh + r*alpha*sinth
ay=-1*r*(omg^2)*sinth + r*alpha*costh
if ax<0.001; %防止除零
ax=0.001;
end;
dir=atan(ay/ax); %求方向
au=cos(dir) + i*sin(dir)
%繪出
linkshape([M(1) M(2)],[x y]); %畫桿
plot([x x+vx],[y y+vy],'color','r','marker','d'); %畫速度
plot([x+vx x+vx+ax],[y+vy y+vy+ay],'color','g','marker','d') %畫加速度
text(x+vx,y+vy,['u=' num2str(nsinth,'%6.1f') '+i' num2str(costh,'%6.1f')]);
text(x+vx+ax,y+vy+ay,['u=' num2str(cos(dir),'%6.1f') '+i' num2str(sin(dir),'%6.1f')]);
%

1.5原點以特定速度,加速度移動

即特定點M以等速水平運動,P速度與加速度方向,和M有加速度時P速度與加速度方向

令原點移動速度為vM=vmx+j*vmy,aM=amx=j*amy
即可利用之前得到的分析結果

Rm = r*cosθ+vmx*t + j*r*sinθ+vmy*t
vm = -r*ω*sinθ+vmx + j*r*ω*cosθ+vmy
am = (-rωcosθ+rαsinθ+amx) +j(-rω2sinθ+rαcosθ+amy)
因此其方向為
uvm=cos(atan((r*ω*cosθ+vmy)/(-r*ω*sinθ+vmx)))+sin(atan((r*ω*cosθ+vmy)/(-r*ω*sinθ+vmx)))
avm= cos(atan((-rωsinθ+rαcosθ)/-rωcosθ+rαsinθ))+j*sin(atan((-rωsinθ+rαcosθ)/-rωcosθ+rαsinθ))


因此我們可以改寫函數如下:
程式架構:
1.資訊
2.求位置
求速度,計算方向
求加速度,計算方向
3.繪製

function mec2_alys(l,theta,omega,alpha,Ori,Oriv,Oria) %分析含有原點移動的參數系統
%作者 b94611001 歐陽太閒
%l=桿長,theta=角度,omega=角速度,alp=角加速度,ori=原點,Oriv=M速度,Oria=M加速度
title('作業10 作者:歐陽太閒 b94611001'); %資訊
xlabel('x-axis cm');
ylabel('y-axis cm');
axis equal;
d2g=pi/180; %角徑轉換
r=l; %桿長
M=Ori; %原點
omg=omega; %角速度
hold on;
theta=theta*d2g;
costh=cos(theta);
sinth=sin(theta);
x=M(1)+ r*costh
y=M(2)+ r*sinth
%求速度
nsinth=-1*sinth;
vx= r*nsinth*omg +Oriv(1)
vy= r*costh*omg +Oriv(2)
dirv=atan(vy/vx); %求方向
vu=cos(dirv) + i*sin(dirv)
%加速度
ax=-1*r*(omg^2)*costh + r*alpha*sinth + Oria(1)
ay=-1*r*(omg^2)*sinth + r*alpha*costh + Oria(2)
if ax<0.001; %防止除零
ax=0.001;
end;
dira=atan(ay/ax); %求方向
au=cos(dira) + i*sin(dira)
%繪出
linkshape([M(1) M(2)],[x y]); %畫桿
plot([x x+vx],[y y+vy],'color','r','marker','d'); %畫速度
plot([x+vx x+vx+ax],[y+vy y+vy+ay],'color','g','marker','d') %畫加速度

例如角度=60,r=5,角速度=1,角加速度=1.2,vM=5+j5,aM=1+j1

mec2_alys(5,60,1,1.2,[0 0],[5 5],[1 1])

d2g =
0.0175
x =
2.5000
y =
4.3301
vx =
0.6699
vy =
7.5000
vu =
0.0890 + 0.9960i
ax =
3.6962
ay =
-0.3301
au =
0.9960 - 0.0890i

結果如圖




由此可知道結果(方向)和不考慮原點M移動有極大不同
以下為疊加vM=0,aM=0的圖形結果


1.6四連桿的運動點P與Q速度與加速度方向與固定桿MR之兩端點之關係
由於本次作業著重於相對速度及加速度的分析,因此同樣利用

[RM]= [AM]= [CM]=
-r3cosθ3 r4cosθ4 ω3 r2ω2cosθ2
-r3sinθ3 r4sinθ4 ω4 r2ω2sinθ2[RM][AM] = [CM]

之關係式求得各桿角速度ω1,ω2,ω3,ω4
注意:此處定義和講義的P,Q定義順序相反,如圖


而我們在前面的分析為M,P,Q,R,以本次分析為準


可求得點P、Q、R等之速度向量，即：

Vp=[r2'+jr2ω2]ejθ2
=jr2ω2(cosθ2+jsinθ2) (1)
Vq=[r2'+jr2ω2]ejθ2+[r3'+jr3ω3]ejθ3
=jr2ω2(cosθ2+jsinθ2)+jr3ω3(cosθ3+jsinθ3) (2)
Ap=[jr2α2-r2ω22]ejθ2
=[jr2α2-r2ω22](cosθ2+jsinθ2) (3)
Aq=[jr2α2-r2ω22]ejθ2 + [jr3α3-r3ω32]ejθ3
=[jr2α2-r2ω22](cosθ2+jsinθ2) + [jr3α3-r3ω32](cosθ3+jsinθ3) (4)

而本次的分析則為(標記上桿序):

桿OP即為1.2,1.3的原點M無速度,加速度桿
Vp=-r*ω*sinθ2 + j*r*ω*cosθ2 (5)
Ap=(-rω2cosθ2+rαsinθ2)+j(-rω2sinθ2+rαcosθ2) (6)
而桿PQ即為1.5的原點M有速度(即p的速度),加速度桿,因此
vm = -r*ω*sinθ+vmx + j*r*ω*cosθ+vmy (7)
Am = (-rωcosθ+rαsinθ+amx) +j(-rω2sinθ+rαcosθ+amy) (8)
將vmx,vmy,vax,vay以vpx,vpy,apx,apy帶入
vq = -r*ω*sinθ3+r*ω*sinθ2 + j*(r*ω*cosθ3+r*ω*cosθ2) (9)
Aq = (-rωcosθ3+rαsinθ3-rω2cosθ2+rαsinθ2) +j(-rω2sinθ3+rαcosθ3-rω2sinθ2+rαcosθ2) (10)


將其展開可以發現有以下關係:
(1)=(5)
(2)=(9)
(3)=(6)
(4)=(10)
而

Vp=-r*ω*sinθ2 + j*r*ω*cosθ2 (5)
Ap=(-rω2cosθ2+rαsinθ2)+j(-rω2sinθ2+rαcosθ2) (6)
vq = -r*ω*sinθ3+r*ω*sinθ2 + j*(r*ω*cosθ3+r*ω*cosθ2) (9)
Aq = (-rωcosθ3+rαsinθ3-rω2cosθ2+rαsinθ2) +j(-rω2sinθ3+rαcosθ3-rω2sinθ2+rαcosθ2) (10)

其方向即為P點對於M(V=0,A=0)的速度與加速度和我們之前所分析的作業8之f4bar.m之分析結果完全相同
故我們可以得知
本次作業和第八次作業有以下共通之處:
1.雖然分析路徑不同,作業八採用的f4bar乃為利用複數法一次將所有值算出,而本次做的則以極座標化為直角座標系的方法(尤拉),將其分析,但

其代數運算轉換後,可以知道兩者是殊途同歸,分析出來的Vp,ap,Vq,aq都是相同的,只是其中q是以相對運動座標系中的旋轉剛體分析出來.
2.其計算,繪圖結果皆同(程式碼計算結果)


pvx =
21.2132
pvy =
21.2132
qvx =
4.1101
qvy =
-24.4963
pax =
-212.1320
pay =
-212.1320
qax =
-1.8712e+003
qay =
-439.0737

二,曲柄滑塊連桿六點瞬心

以下程式利用老師的move4centros為基礎進行修改,可計算在不同曲柄角度時的各桿位置座標
並找尋六點瞬心與曲柄角間之關係.

瞬心之位置須用速度兩點連線之延伸線,可用於求系統上任意點的速度或角速度等,按講義的說法,需按其位置比例求出座標
其實即一個相似三角形求對應位置的程式如下:

因此瞬心分析如下:

function [P,ff]=Xsect(a1,a2,b1,b2)
a1x=real(a1);a1y=imag(a1);a2x=real(a2);a2y=imag(a2);
b1x=real(b1);b1y=imag(b1);b2x=real(b2);b2y=imag(b2);
A=[(a2y-a1y) -(a2x-a1x); (b2y-b1y) -(b2x-b1x)]; %相對座標
determinant=abs(det(A)); %行列式值
if determinant > 10^(-10) %if matrix is not singular %分別處理不同矩陣情況
B=[(a2y-a1y)*a1x-(a2x-a1x)*a1y;...
(b2y-b1y)*b1x-(b2x-b1x)*b1y];
s=A\B; P=complex(s(1),s(2));ff=0; %化為複數型
else % Locate the intersection when the matrix is singular
xx= 10^10;
angle = atan2((a2y-a1y),(a2x-a1x));
P=complex(xx,xx*tan(angle));ff=1;
end


因此可以利用此程式進行瞬心的找尋,
程式結構如下:
1.防呆裝置
2.找邊界
3.建立矩陣和分割點等前置工作
4.代入Xsect
5.找邊界並用sldlink產生資料存入draw4link
6.繪製瞬心軌跡
7.進入for-loop 實際從暫存器中取出資料繪圖

注意:
本程式修改自movef4centrosm,由於其中使用f4bar.m輸出之資料進行運算,但此輸出的矩陣和本次使用的sldlink有所不同,因此需修改輸出矩陣中取得的項位置.

因此測試

>>movesldcentros([0 100 100 0],0,1,0,1,0,5,100)
共可找到六個瞬心,曲桿及滑塊瞬心位於無窮遠.




function movesldcentros(r,th1,td2,tdd2,sigma,driver,ntimes,npts)
%本程式參考movef4centros(r,th1,td2,tdd2,sigma,driver,ntimes,npts)建立
%b94611001 歐陽太閒
%function movef4centros(r,th1,td2,tdd2,sigma,driver,ntimes,npts)
% r: row 滑塊桿資料
% th1: 座標角
% td2,tdd2:角速度
% sigma: 組合模式
% driver: 驅動桿
% ntimes: 迴圈次數
% npts: 分割細緻度

%
%-----------------------------------------------------
clf;
if nargin<10, ntimes=3;npts=100;end; %防呆
figure(1);
[Qstart, Qstop]=sld_angle_limits(r,th1,driver) %先找邊界
npoint=abs(npts);
th2=linspace(Qstart,Qstop,npoint);%分割
val=zeros(6,npoint); %6*npts建立矩陣以存值
for k=1:npoint,
if k>(npoint/2) %克服組合模式sigma=1,=-1的問題
[vr b]=sldlink(r,th1,th2(k),td2,tdd2,sigma,driver); %
else
[vr b]=sldlink(r,th1,th2(k),td2,tdd2,-1*sigma,driver); %
end
% | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
% I | r1 | t1 | td1 | tdd1 | r1' | Vq | |Vq|| ThVq |
% II | r2 | t2 | td2 | tdd2 | r1" | Vp | |Vp|| ThVp |
% III| r3 | t3 | td3 | tdd3 | Q2 | Aq | |Aq|| ThAq |
% IV | r4 | t4 | td4 | tdd4 | P3 | Ap | |Ap|| ThAp |
[pa1,f1]=Xsect(vr(1,1),vr(2,1),vr(4,1),vr(4,1)+10i); %代入Xsect計算瞬心
[pa2,f2]=Xsect(vr(4,1),vr(2,1),vr(1,1),vr(1,1)+10i);
val(1:3,k)=[vr(1,1);vr(2,1);vr(1,1)+vr(4,1)];
val(4:5,k)=[pa2;pa1];
end
%-----------------------------------------------------
x=real(val);y=imag(val); %加入邊界值計算,由於沒有sldlimis,故直接引入sld_angle_limits
[Qstart, Qstop]=sld_angle_limits(r,th1,driver)
[values b]=sldlink(r,th1,Qstart,0,0,1,driver)
if b==1, %可成功組合時
h=draw4link(values,driver); %存入draw4link
line(x(5,:)',y(5,:)','color','r','linestyle','-');
line(x(4,:)',y(4,:)','color','b','linestyle','-.');
range=([1.5*min(min(x(1:3,:))) 2*max(max(x(1:3,:))),...
1.5*min(min(y(1:3,:))) 2*max(max(y(1:3,:)))]);
%axis([-30 30 -30 30]);
axis equal;
grid off;
title('作業10 作者:歐陽太閒 b94611001'); %資訊
%-----------------------------------------------------
for i=2:4,set(h(i),'erasemode','xor','linewidth',6);end %繪圖
h0=line('xdata',[],'ydata',[],'erasemode','xor','color','k','marker','o','linestyle',':');
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',[0 x(5,i) x(3,i) x(4,i) x(1,i)],'ydata',[0 y(5,i) y(3,i) y(4,i) y(1,i)]);
drawnow; %自暫存器讀出
pause(0.1);
end
end % for m loop
%-----------------------------------------------------
end


心得:本次作業第一部份偏重於分析,因此我們有機會能深入的結合機動與動力和微積分的知識剖析多桿系統之運動資訊,而第二部份則利用製圖法求得瞬心,雖然其中為了調整輸出矩陣等問題花了很多時間,但其成果亦相當有趣.

自由軟體界的MATLAB

b94611001 歐陽太閒

Matlab功能強大,擴充模組(toolbox)完整,而且有相當豐富的社群資源(還有一套類似的商業軟體叫做IDL).

可是,如果要在非Windows系統下需要類似的圖形界面套裝軟體時(MATLAB在Mac上因為非原生app,故效率較低),或者在非教育使用等有版權問題的狀況時,就必需找尋同等功能的套裝軟體.

而在自由軟體界,比較常見的有以下兩種,都符合GNU協定,且都有跨平臺支援:

Octave
個人最喜歡的*inx數學運算套件,語法與MATLAB相容,函式庫豐富(可通用MATLAB函式庫),社群資源非常豐富,且也有相當多的toolbox,但要先安裝href>cygwin以模擬GNU environment,使用者多的原因可能是因為Redhat將其作為預設套件

界面:


FreeMat

另一套漂亮好用的自由軟體,不過toolbox較少,但發展的相當快,(短短時間已經到了3.2版)
界面


rlab
原著者Ian Searle的目標是以發展比MATLAB更優美的語法為目標,適合吹毛求疵的power user.
不過可惜的是,目前已經不再更新,但有人接手研發rlabplus


最重要的是,以上軟體都是MATLAB-compatible,也就是相同的語法大多可以通用,就像OO.o那樣不必重新學習.當然,如果需要使用數學符號運算,還有Giac/Xcas,統計可用R等自由軟體套件.

以上是一些個人所知風評不錯的套件,不過茫茫網海,必有遺漏,還請各路高手指正.

Homework 9

Homework 9 b94611001 歐陽太閒
本週(5/3)有來上課,組別A

於題目中可得:
曲桿長Rcm ， 連桿Lcm，滑塊之偏置量為10cm,R=10+1=11,L=R+5=16
無初始速度或加速度,
故可得以下參數:(單位cm)
r(1) =固定桿 =會自行改變設0
r(2) =曲桿 =11
r(3) =結合(連)桿=16
r(4) =垂直偏置量=10
theta1 =第一桿之水平角=0
theta2 = 驅動桿之水平夾角=於計算時代入
td2 = 驅動桿角速度=0
tdd2 = 驅動桿角加速度=0
sigma = 1 組合模數,其中組合模數之閉合或分枝會導致左右鏡像的結果.
driver = 0 驅動桿為第二桿
1 驅動桿為第三桿
2 滑塊驅動


權衡各種分析方法後,使用單純的三角函數處理(如課本)由於反三角函數定義的關系,可能會面臨需撰寫大量if判斷式,且除錯不易的窘境,而Euler化可以大幅簡化計算,就如同電工學的交流電多相合成,複數計算遠較和差化積有效率的多.
故以下作業都將利用老師的sldlink.m(複數計算),sld_angle_limits.m,進行分析
形式：
sldlink(r,theta1,theta2,td2,tdd2,sigma,driver)注解參考附件一
sld_angle_limits(r,theta1,linkdrive) 注解及修改參考附件二
drawsldlimits(r,th1,sigma,driver)
drawsldpaths(r6,th6,r,th1,theta,td2,tdd2,sigma,npts,driver,mode) 注解及修改參考附件三

其中sldlink有部分程式碼在matlab 6.5中運作有問題,即

%f=@(num,ndg);
%round(num*10^ndg)/10^ndg; %取整數
%arg=f(arg,5);

等,但加注解化後對運作無太大影響

一，運動的界限範圍

首先利用老師撰寫的sld_angle_limits代入以下參數進行分析:
1.桿2為驅動桿

r=[0 11 16 10];
theta1=0;
linkdrive=0;
sld_angle_limits(r,theta1,linkdrive);
得到:
Qstart =
3.6000e-006
Qstop =
360.0000
顯然不符,檢查sld_angle_limits程式碼後發現其中缺少了r3≧r4≧0,r3+r4≧r2(即講義7.4.2)之情況
故在其中加入以下程式碼

elseif r3+r4>=r2 & r3>=0 & r3>=r4
Qstart=asin((r4-r3)/r2)
Qstop=pi- asin((r4-r3)/r2)

才能得到:
θmin =θ1 + sin-1((r4-r3) / r2) =asin((10-16)/11)= -0.5769=-33
θmax =θ1 +π- sin-1((r4-r3) / r2) =pi-asin((10-16)/11)=3.7185=213
即
drawsldlimits([0 11 16 10],0,1,0)
Qstart = (起始角)
-33.0557
Qstop = (結束角)
213.0557
其中OP為桿2,PQ為桿3,QR為偏置量
結果如圖

亦可產生

動畫
for k=[-33:213] %簡化程式碼
pause(0.01);
clf;
dwsldshape([0 11 16 10],0,k,1,0)
axis([-16 30 -16 16]); %固定座標系
end



因此:

2.桿3為驅動桿時:
drawsldlimits([0 11 16 10],0,-1,1)
Qstart =
-3.5833
Qstop =
183.5833
結果如圖

3.滑塊為驅動桿時:
drawsldlimits([0 11 16 10],0,-1,2)
Qstart =
-25.0799
Qstop =
25.0799
結果如圖

二，位置分析

由於各點以連桿連結,故參考老師的drawsldlinks.m撰寫如下程式:
其架構如下:
1.資訊
2.將參數放入sldlink且傳回結果
3.複數化
4.若組合存在則繪制
1.接地桿
2.依參數繪桿
3.表示偏心量
4.繪滑塊
5.不能組合時訊息

function [values]=dwsldshape(r,th1,th2,sigma,driver)
%參考老師的drawsldlinks.m
%歐陽太閒 b94611001
%需用到linkshape.m,sldlink.m,sldbox.m
title('作業9.2 作者:歐陽太閒 b94611001');%資訊
xlabel('x-axis');
ylabel('y-axis');
axis equal
grid on
[values b]=sldlink(r,th1,th2,10,0,sigma,driver); %將參數放入sldlink且傳回分析結果
rr=values(:,1);
rr(3)=rr(3)+rr(2);
rx=real(rr);rx(4)=0;%複數化
ry=imag(rr);ry(4)=0;
if b==1
linkshape([0 0],[rx(1) 0]); %繪接地桿
hold on;
linkshape([0 0],[rx(1) ry(1)]); %繪桿2,驅動桿問題已有sldlink處理
linkshape([0 0],[rx(2) ry(2)]); %繪桿3
linkshape([rx(2) ry(2)],[rx(3) ry(3)]);
plot([rx(1) rx(3)],[ry(1) ry(3)],'r-'); % 偏心量,滑塊為重力外力閉合
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('無法組成,角度%6.1f\n',th2); %不能組成時訊息
end


1.桿2為驅動桿(以角度60度為例)

結果如圖


可解析如下
ans =

1.0e+003 *

Columns 1 through 4

桿1位 0.0215 桿1角 0 桿1角速度 0 桿1角加速度 0
桿2位 0.0055 + 0.0095i 桿2角 0.0600 桿2角速度 0.0100 桿2角加速度 0
桿3位 0.0160 + 0.0005i 桿3角 0.0017 桿3角速度 -0.0034 桿3角加速度 0.0599
桿4位 0.0000 + 0.0100i 桿4角 0.0900 桿4角速度 0 桿4角加速度 0

Columns 5 through 8

軸向速度 -0.0936 q速度 -0.0953 + 0.0550i(abs) 0.1100 q速度角度0.1500
軸向加速 -0.7675 p速度 -0.0936 (abs) 0.0936 p速度角度0.1800
Q位置向量 0.0055 + 0.0095i q加速度-0.5500 - 0.9526i(abs) 1.1000 q加速度角度-0.1200
P位置向量 0.0215 + 0.0100i p加速度-0.7675 (abs) 0.7675 p加速度角度 0.1800


2.桿3為驅動桿(以角度60度為例)

結果如圖
dwsldshape([0 11 16 10],0,60,1,1)

ans = (解析同9-二-1)
1.0e+003 *
Columns 1 through 4
0.0183 0 0 0
0.0103 - 0.0039i -0.0205 -0.0078 0.1119
0.0080 + 0.0139i 0.0600 0.0100 0
0.0000 + 0.0100i 0.0900 0 0
Columns 5 through 8
-0.1685 -0.0299 - 0.0800i 0.0854 -0.1105
-0.9896 -0.1685 - 0.0000i 0.1685 -0.1800
0.0103 - 0.0039i -0.1896 + 1.3856i 1.3986 0.0978
0.0183 + 0.0100i -0.9896 0.9896 0.1800

3.滑塊為驅動桿(以角度60度為例)
dwsldshape([0 11 16 10],0,60,1,2)

結果如圖
ans = (解析同9-二-1)
1.0e+002 *
Columns 1 through 4
0 0 0 0
-0.1086 - 0.0175i -1.7085 -0.0100 -0.0108
0.1086 + 0.1175i 0.4725 -0.0100 -0.0016
0.0000 + 0.1000i 0.9000 0 0
Columns 5 through 8
0.1000 -0.0175 + 0.1086i 0.1100 0.9915
0 0.1000 0.1000 0
-0.1086 - 0.0175i 0.0897 + 0.1350i 0.1621 0.5641
0 + 0.1000i 0 0 0

三，路徑分析

而講義7.5提到,可用drawsldpaths.m整合sldlink,body,fb_sld_linits及drawlinks
其中最重要的,便是其可輸出非常漂亮的速度,加速度,角速度,角加速度隨角度變化圖表
以及系統的連續運動路徑

drawsldpaths(r6,th6,r,th1,theta,td2,tdd2,sigma,npts,driver,mode)

r= [0 11 16 10] 各桿長
th1=0 第一桿角度
th2=0 驅動桿或聯結桿角度
td2=0 驅動桿之角速度
tdd2=0 角加速度。
sigma=1 組合模數
driver= 0 曲桿驅動
1 連桿驅動
2 滑塊驅動
Npts=50 分割點數
r6，rh6，nlink：桿上特定點之位置，包括桿長，與桿之夾角及附於何桿。
mode=1 畫所有圖表


由於程式在Matlab6.5下繪制alphapatch會導致畫面混亂,故修改程式,
去除在本分析中用不到的延伸部分,即:

patch('xdata',real(para(1:3)),'ydata',imag(para(1:3)),...
'facecolor','r','facealpha',0.6,'marker','o');hold on;
plot(para(2),'bo');

(請見附件三)

故可得:
1.桿2為驅動桿
drawsldpaths(0,0,[0 11 16 10],0,0.01,0.01,1,50,0,1)

結果如圖
結果如表

2.桿3為驅動桿
drawsldpaths(0,0,[0 11 16 10],0,0.01,0.01,1,50,1,1)

結果如圖
結果如表

3.滑塊為驅動桿
drawsldpaths(0,0,[0 11 16 10],0,0.01,0.01,2,50,2,1)

結果如圖
結果如表

心得:
本次作業為整體性的對一個系統進行分析,故所需用到的程式資源相當的大,且有許多部分需自行改寫,不過也操作的過程中體會MATLAB的強大.

附件一

function [data,form] = sldlink(r,theta1,theta2,td2,tdd2,sigma,driver)

%f=@(num,ndg) round(num*10^ndg)/10^ndg; %取整數,在6.5版會造成無法正確執行的問題
data=zeros(4,8); %生成一個4x8的零矩陣
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 %滑塊驅動桿
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); %避免6.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%桿2或桿3為曲桿
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); %避免6.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, %處理曲桿為3
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)]; %桿2,3的資料交換
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;

附件2

function [Qstart, Qstop]=sld_angle_limits(r,theta1,linkdrive)
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 %分別為講義7.4 RRRP的各種限制
Qstart=asin((r4-r3)/r2);
Qstop=asin((r4+r3)/r2);
elseif r3+r4>=r2 & r4>=r3 %講義程式碼缺少的部分
Qstart=asin((r4-r3)/r2); %r3≧r4≧0,r3+r4≧r2(即講義7.4.2)
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);
elseif r3+r4>=r2 & r3>=0 & r3>=r4
Qstart=asin((r4-r3)/r2);
Qstop=pi- asin((r4-r3)/r2);
else
Qstart=0;Qstop=2*pi;
end
case 1 %coupler %連桿驅動時
if r2-r4<=r3 & r4>=0 & r2>=r4%分別為講義7.4 RRRP的各種限制
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 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 Qstopadjust=(Qstop-Qstart)*1e-8; %避免為零
Qstart=theta1+(Qstart+adjust)*g2d %輸出
Qstop=theta1+(Qstop-adjust)*g2d

附件3

function drawsldpaths(r6,th6,r,th1,td2,tdd2,sigma,npts,driver,mode)

clf;
figure(1);
warning off;
r(abs(r)[Qstart, Qstop]=sld_angle_limits(r,th1,driver) %利用sld_angle_limits取得限制
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); %利用sldlink分析
[para]=body(r6,th6,vr,3); %yongbody畫延伸物(本分析沒有用到)
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 %桿2驅動
val(4:7,i)=[abs(vr(1,1));vr(3,2);vr(3,3);vr(3,4)];
case 1 %桿3驅動
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-');%crank angle or r1曲桿角
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-');%Omega of crank 曲桿角速度
fact=round(max(val(7,:))/max(val(6,:))*10)/10;
hold on;plot(th2,val(7,:)/fact,'b-');%Omega of coupler
xlabel(intitle);ylabel(title3(driver+1));
grid on;
%
subplot(2,2,3);
plot(th2,val(8,:),'k-');% vel of slider 滑塊速度
hold on;plot(th2,val(9,:),'r-');% vel of A
xlabel(intitle);ylabel(title2(1));
grid on;
%
subplot(2,2,4);
plot(th2,val(10,:),'k-');%Acc of slider 滑塊加速度
hold on;plot(th2,val(11,:),'r-');% Acc of A
xlabel(intitle);ylabel(title2(2));
grid on;


第九次作業結束

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
%程式結束