2/4/07

6.4 四連桿分析程式f4bar

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之角加速度為零,因為是輸入值。