四连杆之MATLAB程式讲解.docx

1、四连杆之MATLAB程式讲解第三章 四連桿之MATLAB程式第三章中之四連桿分析可以參考相關資料。本節則針對四連桿之動作程式加以說明。目前所設計之程式有f4bar.m、drawlinks.m、 fb_angle_limits.m、drawlimits.m等四個程式，茲分別說明如下：圖一、四連桿之關係位置及各桿名稱一、f4bar函數： f4bar函數之呼叫格式如下：function values,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)輸入變數：.r(1:4) = 各桿之長度，r(1)為固定桿，其餘分別為曲桿、結合桿及被動桿。 .th

2、eta1 = 第一桿之水平角，或為四連桿之架構角，以角度表示。 .theta2 = 驅動桿之水平夾角，以角度表示。一般為曲桿角，但若為結合桿驅動，則為結合桿之水平夾角。 .td2 = 驅動桿(第二桿或第三桿)之角速度(rad/sec)。 .tdd2 = 驅動桿(第二桿或第三桿)之角加速度(rad/sec2)。 .sigma = +1 or -1. 組合模式，負值表示閉合型，正值為分支型，但有時需視實際情況而定。 .driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) 輸出變數：.form = 組合狀態， 0 ：表示無法組合； 1：可以正確組合 .values = 輸出矩陣，其大小

3、為 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點之速度與加速度量，第五行為向量，第六行為絕對量，第七行為夾角。 值得一提的是第一行、三行、四行及五行之向量表示法屬於複數之型式。故若要得到其絕對值僅需在MATL

4、AB指令檔中，以abs()這一個函數指令即可求得，而以函數angle()則可求得其夾角，雖然第二行與第七行之輸出亦有相對應之夾角。例一：為第二桿為驅動桿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 1 + 1.7321i 2 0 1.8682 + 0.71389i 2 -127.58 173.21

5、 - 100i 200 -236.27 364.19 - 953.09i 1020.3 Column 7 60 20.913 -30 -69.087 form = 1 (表示可以組合)本例中，有框線者表示其為輸入值。但第一行則已經轉換為複數型式。未來複數型式要轉為x-y座表時，只要使用函數real()及imag()兩指令，即可進行轉換。例二：為第三桿(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.46

6、41i 60 10 0.33205 + 1.9722i 80.443 24.333 Columns 4 through 6 0 1.3321 - 1.4919i 2 -582.55 0.33205 + 1.9722i 2 0 -988.55 - 882.66i 1325.3 496.46 188.64 - 31.759i 191.29 Column 7 -48.239 80.443 -138.24 -9.5568 form = 1程式內容：function values,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)%function v

7、alues,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver)% program designed by Din-sue Fon, NTU, revised from Waldrons% This function analyzes a four-bar linkage when the crank is the % driving link. The input values are:% theta1,theta2 are angles in degrees%r(1) = length of vector 1 (frame)%r(2) =

8、length of vector 2 (crank)%r(3) = length of vector 3 (coupler)%r(4) = length of vector 4 (rocker or slider offset)%td2 = crank or coupler angular velocity (rad/sec)%tdd2 = crank or coupler angular acceleration (rad/sec2)%sigma = +1 or -1. Identifies assembly mode%driver = 0 for crank as driver; 1 fo

9、r coupler as driver% The results are returned in the vector values. The answers are % stored in values according to the following:%values (1:4,1) = link position%values (1:4,2) = link angles in degrees%values (1:4,3) = link angular velocities%values (1:4,4) = link angular accelerations%values (1,5)

10、= velocity of point Q%values (2,5) = velocity of point P%values (3,5) = acceleration of point Q%values (4,5) = acceleration of point P%vakyes (4,6) =absolute values of values(:,4)%vakyes (4,7) =angles in degrees of values(:,4)%form = assembly flag. If form = 0, mechanism cannot be% assembled.%conver

11、t input datavalues=zeros(4,7);% if coupler is the driver, interchange the vetor 3 & 2% if theta2180|theta2=0) form=1;% Check for the denominator equal to zero if abs(C-A)=eps t4=2*atan(-B+sigma*sqrt(arg)/(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=s

12、in(t3); c3=cos(t3); elseif abs(C-A)=0) t3=2*atan(-B-sigma*sqrt(arg)/(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

13、)*td(2)*sx;r(2)*td(2)*cx; CM=AMBM; 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=AMBM; tdd(3)=CM(1); tdd(4)=CM(2); %store

14、results in array values % coordinates of P and Q if driver=1, r=r(1) r(3) r(2) r(4); c2=c3;c3=cx;s2=s3;s3=sx; td=td(1) td(3) td(2) td(4); tdd=tdd(1) tdd(3) tdd(2) tdd(4); theta=theta(1) theta(3) theta(2) theta(4); else c2=cx;s2=sx; end for j=1:4, values(j,1)=r(j).*exp(i*theta(j); values(j,2)=theta(j

15、)/fact; values(j,3)=td(j); values(j,4)=tdd(j); end % position vectors values(1,5)=r(2).*exp(i*theta(2);%velocity for point Q values(2,5)=r(4).*exp(i*theta(4);%velocity for point P values(3,5)=i*r(2).*(tdd(2)-td(2).*td(2).*exp(i*theta(2);%accel of Q values(4,5)=i*r(4).*(tdd(4)-td(4).*td(4).*exp(i*the

16、ta(4);%accel of P for j=1:4, values(j,6)=abs(values(j,5); %absolute values for values(:,4) values(j,7)=angle(values(j,5)/fact; %angles for values(:,4) end %find the accelerationselse form=0; if driver=1, r=r(1) r(3) r(2) r(4); for j=1:4, values(j,1)=r(j).*exp(i*theta(j);end % position vectors endend

17、二、drawlinks函數 drawlinks之目的在利用MATLAB繪製四連桿之相關位置，其本身會呼叫f4bar.m函數以計算四連桿之向量位置，然後繪圖。其呼叫格式如下： function drawlinks(r,th1,th2,sigma,driver) 其輸入各式與f4bar.m大體相同，茲說明如下：.r(1:4) = 各桿之長度，r(1)為固定桿，其餘分別為曲桿、結合桿及被動桿。 .theta1 = 第一桿之水平角，或為四連桿之架構角，以角度表示。 .theta2 = 驅動桿之水平夾角，以角度表示。一般為曲桿角，但若為結合桿驅動，則為結合桿之水平夾角。 .sigma = +1 or -

18、1. 組合模式，負值表示閉合型，正值為分支型，但有時需視實際情況而定。 .driver = 0 (驅動桿為第二桿); 1 (驅動桿為第三桿) 例三、第二桿為驅動桿 drawlinks(3 2 4 2,0,60,-1,0)圖二、四連桿之繪圖其繪出之四連桿為如圖二。黑色為第一桿，藍色為第二桿，紅色為第三桿，綠色為第四桿。例二、第三桿為結合桿(coupler) drawlinks(3 2 4 2,0,60,-1,1)圖三、以結合桿為驅動桿(r =3 2 4 2) 圖三即為所得之答案，此時四連桿為分支型(branch)，因為目前之情況無法轉為閉合型，即使將sigma值變號，仍為分支型，如： clf；d

19、rawlinks(3 2 4 2,0,60,1,1)圖四、當sigma=1時並以結合桿為驅動桿(r =3 2 4 2) 利用drawlinks亦可繪出各種角度之圖型，可以作為四連桿運動過程之觀察，相當方便，例如：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,0);end圖五、多重位置之四連桿運動情形(r =1 3 3 4)其結果如圖五，當然若以第三桿為驅動桿時，亦可獲得同樣的結果，例如： clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,1);end圖六、以第三桿為驅動桿之情形(r =1 3 3 4)當sigm

20、a 變號時，亦可看出其不同的轉動方式，如下：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,1,1);end圖七、當sigma變號(為正時) (r =1 3 3 4)drawlinks程式內容function drawlinks(r,th1,th2,sigma,driver)%function drawlinks(r,th1,th2,sigma,driver)%draw the positions of four-bar links %will call f4bar.m funcion%designed by Din-Sue Fon, NTU%r: row

21、vector for four links%th1: frame angle%th2: crank angle or couple angle%sigma: assembly mode%driver: 0 for crank, 1 for coupler%clf;r b=f4bar(r,th1,th2,0,0,sigma,driver);r(3,1)=r(1,1)+r(4,1);rx=real(r(:,1);rx(4)=0;ry=imag(r(:,1);ry(4)=0;if b=1 plot(0 rx(1),0 ry(1),k-,LineWidth,4); hold on; if driver

22、=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); plot(rx,ry,bo); text(0,0, O);text(rx(1),ry(1), R); text(rx(2),ry(2), P);text(rx(3)

23、,ry(3), Q);else fprintf(Combination of links fail at degrees %6.1fn,th2);endaxis equalgrid on三、drawlimits函數四連桿之迴轉過程，能完全迴轉的情況仍然很少，有些時候無法獲得完整一圈的迴轉。亦即依葛列夫定理四連桿之第一或第二類類型決定，前者為完整迴轉型，後者則有迴轉角度之限制，這些限制因四連桿長度決定之。四連桿迴轉過程中，有可能其中兩桿會連成一線，或重疊成一線，前者若成立時，即變成三角形，後者若重疊時，亦會構成另一個三角形。理論上連桿構成三角形應不會有相對運動。故可稱為四連桿之運動極限。由這兩個

24、極端位置，可以知道四連桿之最終運動限制。A. 第二桿為驅動桿在數學上，表示這兩個狀況之方法可以利用下列二種不等式進行測試：r1+r2|r3-r4|而由其不等式之方向，可以構成四種狀況，並進而求得該狀況之角度。下面為第二桿為驅動桿時之四種情況：(1) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min =0 ，max =2(2) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min = - cos-1r12+r22-(r3+r4)2/(2 r1 r2)，max = cos-1r12+r22-(r3+r4)2/(2 r1 r2)(3) 當r1+r2r3+r4，|r1-r2|r3-

25、r4| 時min = cos-1r12+r22-(r3-r4)2/(2 r1 r2)，max = cos-1r12+r22-(r3+r4)2/(2 r1 r2)(4) 當r1+r2r3+r4，|r1-r2|r3-r4| 時min = cos-1r12+r22-(r3-r4)2/(2 r1 r2) ，max = 2- cos-1r12+r22-(r3-r4)2/(2 r1 r2)B. 第三桿為驅動桿第三桿結合桿為驅動桿時，則仍然取決於四連桿屬於葛列斯荷(Grashof)一型或二型。若屬一型連桿，則當第三桿r3為最短桿時，第三桿可以作360度迴轉。其餘之限制條件雖不如以第二桿為驅動桿者，但其極;限

26、狀況是當第二桿與第四桿相平行時，變成無法繼續迴轉，除非它是處於平行四邊形。將四種情況依下列二不等式之情況加以分類，在這些分類中，若兩式均為等號時，則應歸屬於第五類：r1+r3|r2-r4|(5) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min =0 ，max =2(6) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min = - cos-1r12+r32-(r2+r4)2/(2 r1 r3) ，max = cos-1r12+r32-(r2+r4)2/(2 r1 r3)(7) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min = cos-1r12+r32-(r

27、2-r4)2/(2 r1 r3) ，max = cos-1r12+r32-(r2+r4)2/(2 r1 r3)(8) 當r1+r3r2+r4，|r1-r3|r2-r4| 時min = cos-1r12+r32-(r2-r4)2/(2 r1 r3) ，max = 2- cos-1r12+r32-(r2-r4)2/(2 r1 r3)C. fb_angle_limits函數觀察上面討論之四個極限角度，可以寫一組程式進行計算。由於以第三桿驅動與第二桿驅動，在計算上僅是將其中之r2與r3之位置對調即可。為尋找上述極限角度min、max，可用函數fb_angle_limits進行尋找，其格式如下：function Qstart, Qstop=fb_angle_limits(r,Q1,driver) 其中輸入項目有： r = 四連桿之長度向量，其定義與前函數相同。 Q1 = 第一桿之夾角，角度表示(deg)。 driver = 驅動模式(=0 第二桿驅動； =1 第三桿驅動)。而輸出項為兩個角度：Qstart = 驅動桿(第二桿或第三桿)之最低限角度 (d