2007年6月16日 星期六

作業八

B94611003

一.
設桿2角度theta2=45度時,求各點之位置、速度與加速度為何?

本題採用講義上之f4bar函式
分析出四連趕各桿位置,速度,加速度,然後取值

以下為函式的內容:

function [data,form] = f4bar(r,theta1,theta2,td2,tdd2,mode,linkdrive)
%輸入參數為
%r= [r1 r2 r3 r4]: 各桿之長度
%theta1= 桿一之水平角度
%theta2= 驅動桿(可為桿2或3)之水平角度
%td2= 驅動桿之角速度(rad/sec)
%tdd2= 驅動桿之角加速度(rad/sec^2)
%mode= +1 or -1 組合模數,負值表閉合型;正值表分支型
%linkdrive = 0表示驅動桿為第二桿; 1表示驅動桿為第三桿
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([4 3 3 5],0,45,10,0,-1,0)

data =

1.0e+003 *

Columns 1 through 5

0.0040 0 0 0 0.0212 + 0.0212i
0.0021 + 0.0021i 0.0450 0.0100 0 0.0041 - 0.0245i
0.0011 + 0.0028i 0.0695 -0.0163 0.4914 -0.2121 - 0.2121i
-0.0008 + 0.0049i 0.0995 -0.0050 0.3836 -1.8712 - 0.4391i

Column 6

0.0021 + 0.0021i
0.0032 + 0.0049i
0
0

function dyad_draw([4 3 3 5],[0 45 69 99],[0 10 16.3 5],[0 0 491.4 383.6])

%column1: 各桿之位置向量(複數型式)
%column2: 各桿之水平夾角
%column3: 各桿之角速度
%column4: 各桿之角加速度
%column5-1: Q點之速度
%column5-2: P點之速度
%column5-3: Q點之加速度
%column5-4: P點之加速度
%column6-1: Q點之位置向量
%column6-2: P點之位置向量

form =

1

%form=1 表示此四連桿可構成,若為0 則表不行


由程式結果得到
O點位置: 0 + 0i
Q點位置: 0.0021 + 0.0021i
P點位置: 0.0032 + 0.0049i
R點位置: 0.0040 + 0i(同桿一)


>> abs(data(:,5))

ans =

1.0e+003 *

0.0300
0.0248
0.3000
1.9220

又對column5取絕對值後
Q點速度為0.03
P點速度為0.0248
R點和O點速度為0(因桿1為固定桿)
(速度單位為: 長度/時間)

同理
Q點加速度為0.3
p點加速度為1.9220
R點和O點加速度為0
(加速度為: 長度/時間^2)

第一桿:0.0040
第二桿:0.0021 + 0.0021i
第三桿:0.0011 + 0.0028i
第四桿:-0.0008 + 0.0049i
速度如下:
第一桿:0(rad/s)
第二桿:10.0000(rad/s)
第三桿:16.2681(rad/s)
第四桿:4.9677(rad/s)
加速度:
第一桿:0(rad/s*s)
第二桿:0(rad/s*s)
第三桿:491.4428(rad/s*s)
第四桿:383.6120(rad/s*s)

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

本題使用講義上function drawlinks函式
繪出給定條件drawlinks([4 3 3 5],0,45,1,0)之四連桿

以下為函式內容:
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 mode="1;end" rr="values(:,1);" rx="real(rr(:,1));rx(4)=" ry="imag(rr(:,1));ry(4)=" b="=" linkdrive="=">> axis([-1 5 -1 5]);
axis equal;
drawlinks([4 3 3 5],0,45,-1,0)
%在顯示畫面上標明三點速度方向及大小
gtext ('速度方向0.0212 + 0.0212i 大小0.03)
gtext ('速度方向0.0041 - 0.0245i 大小0.0248)
gtext ('速度方向0 + 0 大小0')
gtext ('速度方向0 + 0 大小0')


ans =

Columns 1 through 5

4.0000 0 0 0 0
2.1213 + 2.1213i 45.0000 0 0 0
1.0513 + 2.8098i 69.4856 0 0 0
-0.8274 + 4.9311i 99.5246 0 0 0

Column 6

2.1213 + 2.1213i
3.1726 + 4.9311i
0
0

%Q速度方向0.0212 + 0.0212i 大小0.03
%P速度方向0.0041 - 0.0245i 大小0.0248
%R及0速度方向0 + 0 大小0

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

老師第六章講義function drawlimits在角度限制上的運用下
求出當桿2迴轉時,此組四連桿之限制角度
再運用fb_angle_limits求之

function [theta]=deadangle(r,mode)
% Finding the max. and min. angles of link 4.
% Input: r: linkage matrix.
% Output:theta:angles of link 2 & link 4.
% Examples:
rr=r.*r;d2g=pi/180;if nargin<2,mode=1;end ph1="acos(((r(2)-r(3))^2-rr(1)-rr(4))/(2*r(1)*r(4)));" ph2="acos(((r(2)+r(3))^2-rr(1)-rr(4))/(2*r(1)*r(4)));" th1="acos((rr(1)+rr(2)-(r(3)+r(4))^2)/(2*r(1)*r(2)));" th2="acos((rr(1)+rr(2)-(r(3)-r(4))^2)/(2*r(1)*r(2)));" th1="0;end" th2="2*pi;end">th2,temp=th1;th1=th2;th2=temp;end
if th1==pi&th2==2*pi,th1=0;end
if th1==0&th2==pi,th2=2*pi;end
if th1==0&th2==0,th2=2*pi;end
if ~isreal(ph1),ph1=0;end
if ~isreal(ph2),ph2=2*pi;end
if ph1>ph2,temp=ph1;ph1=ph2;ph2=temp;end
if ph1==0&ph2==0,ph2=2*pi;end
theta=[th1,th2;ph1,ph2]/d2g;


將theta=deadangle([4 3 3 5])執行後
可以得到:
theta =

0 28.9550
0 97.1808

上死點:















四.

這三個角度是無法繪製出來的,因theta2只能介於28.9550與97.1808之間,

三個角度都超過死點,所以東西出不來。















五.

需要用到下列函式(drawlinkscg.m)

function [values]=drawlinks(r,th1,th2,mode,linkdrive)
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
title('Hit Ctrl+C to stop');
axis([-5 5 -5 5]);
grid on

以主程式來呼叫
閉合型
for t=1:6
for i=29:331
clf;drawlinkscg([4 3 3 5],0,i,-1,0);
pause(0.0000001);
end;for i=29:331
clf;drawlinkscg([4 3 3 5],0,360-i,-1,0);
pause(0.0000001);
end;

分支型
for t=1:6
for i=29:331
clf;drawlinkscg([4 3 3 5],0,i,1,0);
pause(0.0000001);
end;for i=29:331
clf;drawlinkscg([4 3 3 5],0,360-i,1,0);
pause(0.0000001);
end;

沒有留言: