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迴圈所製成的動畫

2 則留言:

Martin Foon 提到...

做得很好,最後一個的範圍應設得大一些,才能將運動的內容穩住。且 axis之設限指令應拿到迴圈之外。

大鳥 提到...

我跟老師一樣,座標軸的範圍可以找適當一點,這樣圖比較漂亮!!