作業八(改)

8A: 本人本週四(4/26)曾來上課
有一組四連桿，其桿長分別為r=[4 3 3 5]，由桿2驅動，設第一固定桿角度theta1=0度; 角速度
td2=10rad/s; 角加速度tdd2=0 rad/s^2。

8B.1:
設桿2角度theta2=45度時，求各點之位置、速度與加速度為何？

先建立一程式
f4bar.m:
function [data,form]=f4bar(r,theta1,theta2,v2,a2,mode,linkdrive)
% This function analyzes a four-bar linkage when the driving link is
% crank(r2曲柄啟動) or coupler(r3聯結桿啟動). The input data are:theta1,theta2 are angles in degrees
% r=[r1 r2 r3 r4]= lengths of links(1=frame)
%v2 = crank or coupler angular velocity (rad/sec)
%a2 = 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(r2端點)
%data (2,5)= velocity of point P(r3端點)
%data (3,5)= acceleration of point Q(r2端點)
%data (4,5)= acceleration of point P(r3端點)
%data (1,6)= position of Q(r2端點)
%data (2,6)= position of P(r3端點)
%form = assembly status. form = 0, mechanism fails to assemble.
if nargin<7,linkdrive=0;end
if nargin<6,mode=1;end
data=zeros(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
rr=r.*r;d2g=pi/180;
[theta,v,a]=deal(zeros(4,1));
theta(1:2)=[theta1 theta2]*d2g;
t1=theta(1);tx=theta(2);
s1=sin(t1);c1=cos(t1);
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,
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
v(2)=v2;
AM=[-r(3)*s3, r(4)*s4; -r(3)*c3, r(4)*c4];
BM=[r(2)*v(2)*sx;r(2)*v(2)*cx];
CM=AM\BM;
v(3)=CM(1);v(4)=CM(2);
%acceleration calculation
a(2)=a2;
BM=[r(2)*a(2)*sx+r(2)*v(2)*v(2)*cx+r(3)*v(3)*v(3)*c3-...
r(4)*v(4)*v(4)*c4;r(2)*a(2)*cx-r(2)*v(2)*v(2)*sx-...
r(3)*v(3)*v(3)*s3+r(4)*v(4)*v(4)*s4];
CM=AM\BM;
a(3)=CM(1);a(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)];
v(2:3)=[td(3) v(2)];a(2:3)=[a(3) a(2)];
else
c2=cx;s2=sx;
end
for j=1:4,
data(j,1:4)=[r(j)*exp(i*theta(j)) theta(j)/d2g v(j) a(j)] ;
end % position vectors
data(1,5)=r(2)*v(2)*exp(i*theta(2));%velocity for point Q
data(2,5)=r(4)*v(4)*exp(i*theta(4));%velocity for point P
data(3,5)=r(2)*(i*a(2)-v(2)*v(2))*exp(i*theta(2));%acc of Q
data(4,5)=r(4)*(i*a(4)-v(4)*v(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
接著執行
[val,form]=f4bar([4 3 3 5],0,45,10,0,-1,0)
先求四桿之角度
abs(val(:,2))'
ans =
0 45.0000 69.4856 99.5246
以r1-r2接點為O,r2-r3:P,r3-r4:Q,r4-r1:R
r2x=3*cosd(45)
r2x=
2.1213
r2y=3*sind(45)
r2y=
2.1213
r3x=3*cosd(45)+3*cosd(69.4856)
r3x =
3.1726
r3y=3*sind(45)+3*sind(69.4856)
r3y=
4.9311

各端點位置:
O(0,0)
R(4,0)
P(r2x,r2y)=(2.1213,2.1213)
Q(r3x,r3y)=(3.1726,4.9311)
-------------------------------Ans

先求四桿之角速度
abs(val(:,3))'
ans =
0 10.0000 16.2681 4.9677

各端點速度:
VO=0
VP=3*10=30
VQ=5*4.9677=24.8385
VR=0
----------------------------Ans

先求四桿之角加速度
abs(val(:,4))'
ans =
0 0 491.4428 383.6120
ar=r''+r*(theta)^2
atheta=r*theta''+2*r'*theta'
A=(ar^2+atheta^2)^0.5
Par=3*(10)^2=300
Patheta=0
AP=300
Qar=5*(4.9677)^2=123.3902
Qatheta=5*383.6120=1.9181e+003
AQ=((123.3902)^2+(5*383.6120)^2)^0.5=1.9220e+003

各端點加速度:
AO=0
AP=300
AQ=1.9220e+003
AR=0
--------------------------Ans

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

先建立一程式:f4bar.m
再執行[val,form]=f4bar([4 3 3 5],0,45,10,0,-1,0);
q=1:4;
Vx(q)=abs(val(q,1)).*abs(val(q,3)).*sind(abs(val(q,2)));
Vy(q)=abs(val(q,1)).*abs(val(q,3)).*cosd(abs(val(q,2)));
求出
Vx(2)=21.2132
Vy(2)=21.2132
Vx(4)=24.4963
Vy(4)=-4.1101
再執行
r=[4 3 3 5];
th1=0;
th2=45;
mode=-1;
linkdrive=0;
r=[r(1) r(3) r(2) r(4)];
[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,V=(0,0)');text(rx(1),ry(1),' R,V=(0,0)');
text(rx(2),ry(2),' P,V=(21.2132,21.2132)');text(rx(3),ry(3),'Q,V=(24.4963,-4.1101)');
else
fprintf('Combination of links fail at degrees %6.1f\n',th2);
end
axis equal
grid on
圖示


8B.3
當桿２迴轉時，求出此組四連桿之限制角度，並繪出其位置（以程式為之）。

先建ㄧ程式:deadangle.m
deadangle.m:
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)));
if ~isreal(th1),th1=0;end
if ~isreal(th2),th2=2*pi;end
if th1>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
得r2之限制角為28.9550
r4之限制角為97.1808
以閉合型及分支型區分
閉合型:

分支型:

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

閉合型:

分支型:

有些沒有值是因為theta2已超過死點之角度,請看8B.3
8B.5
若將問題三考慮在內，只在可迴轉的範圍內迴轉，請問你能讓此組四連桿作成動畫方式迴轉嗎？

先建立一程式:drawlinkscg.m
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:5
for i=29:331
clf;drawlinkscg([4 3 3 5],0,i,-1,0);
pause(0.00000000001);
end;for i=29:331
clf;drawlinkscg([4 3 3 5],0,360-i,-1,0);
pause(0.00000000001);
end;
end
動畫

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