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


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