6.4程式介紹
程式參數
以MATLAB寫作，可以獲得以上所述之解法，下面為f4bar()函數之內容，係根據上面的分析過程撰寫而成。其呼叫之型式如下，其中輸入有七項，分別說明如下：
[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
其中第一行之連桿位置向量，屬於單桿的位置向量，以格式以複數表示。第二行為各桿之水平夾角，以度表示；第三及第四行為各桿之角度速度及角加速度，以單位時間之弧度表示。第五至七行則為P點與Q點之速度與加速度量，第五行為向量，第六行為絕對量，第七行為夾角，以度數表示。
值得一提的是第一行之向量表示法屬於複數之型式，故若要得到其絕對值僅需在MATLAB指令檔中，以abs()這一個函數指令即可求得，而以函數angle()則可求得其夾角，雖然第二行與第七行之輸出亦有相對應之夾角。
其中第一行之連桿位置是比較爭議性的問題，目前有不同版本，有些是單桿的位置向量。但第五至七行之速度與加速度量則為絕對量值。
程式內容
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);
% 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,td,tdd]=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
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
執行範例
(1)桿2驅動：
桿2驅動，r=[3 2 4 2]，固定桿角度theta1=0度; 桿2角度theta2=60度;角速度 td2=10rad/s; 角加速度tdd2=0rad/s^2; 組合模數mode=-1; driver=0(crank)
>> [val,form]=f4bar([3 2 4 2],0,60,10,0,-1,0)
val =
Columns 1 through 3
3 0 0
1 + 1.7321i 60 10
3.8682 - 1.0182i -14.746 5.4078
1.8682 + 0.71389i 20.913 16.549
Columns 4 through 6
0 10 + 17.321i 1 + 1.7321i
0 30.919 + 11.815i 4.8682 + 0.71389i
-127.58 -100 - 173.21i 0
-236.27 -343.02 - 636.93i 0
form =
1
上述之結果，form = 1表示可以構成四連桿。特別注意第一行資料雖為位置向量，但仍以複數型式表示，以期能得到每一桿之向量值，以供其他程式呼叫使用。其絕對值仍應為[3 2 4 2]。如：
>> abs(val(:,1))'
ans = 3 2 4 2
四連桿之角度為：
>> abs(val(:,2))'
ans = 0 60 14.746 20.913
桿1之角度為零，桿2為60度，符合已知之資料。
各桿之速度則為
>> abs(val(:,3))'
ans = 0 10 5.4078 16.549
同樣桿1之角速度為零，而桿2為10rad/s，與原已知相同。
各桿之加速度為：
>> abs(val(:,4))'
ans = 0 0 127.58 236.27
桿1及桿2之加速度為零，而桿3為127.58 rad/s^2，桿4為236.27rad/s^2。
(1)桿3驅動：
第三桿(coupler)為驅動桿，r=[3 2 4 2]，固定桿角度theta1=0度; 桿3角度theta2=30度; 桿3角速度td2=10rad/s;桿3角加速度 tdd2=0rad/s^2;模數 mode=-1; 驅動桿driver=1(coupler)。
>> [val,form]=f4bar([3 2 4 2],0,60,10,0,-1,1)
val =
Columns 1 through 3
3 0 0
1.3321 - 1.4919i -48.239 -8.9487
2 + 3.4641i 60 10
0.33205 + 1.9722i 80.443 24.333
Columns 4 through 6
0 -11.92 + 13.35i 1.3321 - 1.4919i
-582.55 8.0799 + 47.991i 3.3321 + 1.9722i
0 -975.76 - 656.52i 0
496.46 -1175.8 - 1002.9i 0
form =
1
同樣，各桿之值為：
>> abs(val(:,1))'
ans =
3 2 4 2
各桿相對應之角度如下：
>> abs(val(:,2))'
ans = 0 48.239 60 80.443
其中桿3之角度為60度，因為該桿為輸入值。
四連桿之角速度分別為：
>> abs(val(:,3))'
ans = 0 8.9487 10 24.333
桿3之角速度為10rad/s，與輸入相同，此時桿2之角速度則為8.9487 rad/s。
各桿之角加速度分別為：
>> abs(val(:,4))'
ans = 0 582.55 0 496.46
桿3之角加速度為零，因為是輸入值。
No comments:
Post a Comment