第三章四連桿之MATLAB程式

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

2、角，或為四連桿之架構(gòu)角，以角度表示。.theta2 =驅(qū)動桿之水平夾角，以角度表示。一般為曲桿角，但若為結(jié)合 桿驅(qū)動，則為結(jié)合桿之水平夾角。.td2 =驅(qū)動桿(第二桿或第三桿)之角速度(rad/sec)。.tdd2 =驅(qū)動桿(第二桿或第三桿)之角加速度(rad/secA2)o.sigma =+1 or -1.組合模式，負值表示閉合型，正值為分支型，但有時需視實際情況而定。.driver0 (驅(qū)動桿為第二桿)；1 (驅(qū)動桿為第三桿)輸出變數(shù):.form組合狀態(tài)，0 :表示無法組合；1:可以正確組合.values =輸出矩陣，其大小為4 X 7 ，各行之資料分配如下:12 (deg)3(rad/s

3、)4(rad/s 2)567I桿1位置9 13 1a 1VQ|VQ|/ VQ :II桿2位置9 23 2a 2VP|VP|/ VPIII桿3位置9 33 3a 3AQ :|AQ|/ AQ :IV桿4位置9 43 4a 4AP|AP|/ AP其中第一行之連桿位置向量，屬於單桿的位置向量。第二行為各桿之水平夾 角，第三及第四行為各桿之角度速度及角加速度。第五至七行則為P點與Q點之速度與加速度量，第五行為向量，第六行為絕對量，第七行為夾角。值得一提的是第一行、三行、四行及五行之向量表示法屬於複數(shù)之型式。故 若要得到其絕對值僅需在 MATLAB指令檔中，以abs()這一個函數(shù)指令即可求 得，而以函數(shù)a

4、ngle()則可求得其夾角，雖然第二行與第七行之輸出亦有相對 應(yīng)之夾角。例一：為第二桿為驅(qū)動桿val,form=f4bar(3 2 4 2,0,60,10,0,-1,0)val = Columns 1 through 331 +1.7321i3.8682 -1.0182i1.8682 +0.71389iColum ns 4 through 600-127.58-236.27Colu mn 76020.913-30-69.087form =1 (表示可以組合)060-14.74620.9130105.407816.5491 +1.7321i21.8682 +0.71389i2173.21 -10

5、0i200364.19 -953.09i1020.3本例中，有框線者表示其為輸入值。但第一行則已經(jīng)轉(zhuǎn)換為複數(shù)型式。未來複數(shù) 型式要轉(zhuǎn)為x-y座表時，只要使用函數(shù)real()及imag()兩指令，即可進行轉(zhuǎn)換。例二：為第三桿（coupler） 為驅(qū)動桿val,form=f4bar(3 2 4 2,0,60,10,0,-1,1)val =Columns 1 through 331.3321 -2 +0.33205 +1.4919i3.4641i1.9722i0 -48.239 60 80.443-8.94871024.333Columns 4 through 60-582.550496.461.3

6、321 -1.4919i20.33205 +1.9722i-988.55 -882.66i1325.3188.64 -31.759i191.29Colu mn 7-48.23980.443-138.24-9.5568form =1程式內(nèi)容:function values,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver) %function values,form = f4bar(r,theta1,theta2,td2,tdd2,sigma,driver) % program designed by Din-sue Fon, NTU, rev

7、ised from Waldron's % 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)%r(2)%r(3)%r(4)%td2 %tdd2 %sigma=length of vector 1 (frame)=length of vector 2 (crank)=length of vector 3 (coupler)=length of vecto

8、r 4 (rocker or slider offset)=crank or coupler angular velocity (rad/sec)=crank or coupler angular acceleration (rad/secA2)=+1 or -1. Identifies assembly mode%driver = 0 for crank as driver; 1 for coupler as driver% The results are returned in the vector "values". The answers are % stored

9、in values according to the following:%values (1:4,1)%values (1:4,2)%values (1:4,3)%values (1:4,4)=link position=link angles in degrees=link angular velocities=link angular accelerations%values (1,5) %values (2,5) %values (3,5) %values (4,5) %vakyes (4,6) %vakyes (4,7) %form%= velocity of point Q= ve

10、locity of point P= acceleration of point Q= acceleration of point P =absolute values of values(:,4) =angles in degrees of values(:,4)= assembly flag. If form = 0, mechanism cannot be assembled.%convert input datavalues=zeros(4,7);% if coupler is the driver, interchange the vetor 3 & 2% if theta2

11、>180|theta2<0, sigma=-sigma;endif driver=1,r=r(1) r(3) r(2) r(4);endrr=r.*r; fact=pi/180; theta=zeros(4,1);td=zeros(4,1); tdd=zeros(4,1);theta(1:2)=theta1 theta2*fact; 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

12、=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; arg=B*B-C*C+A*A;if (arg>=0)form=1;% Check for the denominator equal to zeroif 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)

13、; s3=sin(t3);c3=cos(t3);elseif abs(C-A)<eps% If the denominator is zero, compute theta(3) firstA=-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); arg=B*B-C*C+A*A;if (arg>=0) t3=2*atan(-B-sigma*sqrt(arg)/(C-A); s3=sin(t3); c3=

14、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);endendtheta(3)=t3;theta(4)=t4;%velocity calculationtd(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=AMBM;td(3)=CM(1);td(4)=CM(2);%acceleration calculationtdd(2)=tdd2; BM=r

15、(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 results in array values% coordinates of P and Qif driver=1,r=r(1) r(3) r(2) r(4); c2=c3;c3=cx;s2=s3;s3=sx; t

16、d=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);elsec2=cx;s2=sx;endfo門=1:4,values(j,1)=r(j).*exp(i*theta(j); values(j,2)=theta(j)/fact; values(j,3)=td(j);values(j,4)=tdd(j);end % position vectorsvalues(1,5)=r(2).*exp(i*theta(2); %velocity for poi

17、nt 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 Qvalues(4,5)=i*r(4).*(tdd(4)-td(4).*td(4).*exp(i*theta(4);%accel of Pfor j=1:4,values(j,6)=abs(values(j,5); %absolute values for values(:,4) values(j,7)=angle(values(j,5)/f

18、act; %angles for values(:,4) end%find the accelerationselseform=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 vectorsendend1、drawlinks 函數(shù)drawlinks之目的在利用MATLAB繪製四連桿之相關(guān)位置，其本身會呼叫f4bar.m函數(shù)以計算四連桿之向量位置，然後繪圖。其呼叫格式如下:function drawlinks(r,th1,th2,sigma,driver

19、)其輸入各式與f4bar.m大體相同，茲說明如下:.r(1:4)=各桿之長度，r(1)為固定桿，其餘分別為曲桿、結(jié)合桿及被動桿。.thetal =第一桿之水平角，或為四連桿之架構(gòu)角，以角度表示。.theta2 =驅(qū)動桿之水平夾角，以角度表示。一般為曲桿角，但若為結(jié)合 桿驅(qū)動，則為結(jié)合桿之水平夾角。.sigma =+1 or -1.組合模式，負值表示閉合型，正值為分支型，但有時需視實際情況而定。.driver = 0 (驅(qū)動桿為第二桿);1 (驅(qū)動桿為第三桿)例三、第二桿為驅(qū)動桿drawlinks(3 2 4 2,0,60,-1,0)2.51.50.5-0.51.52.53.54.57f/瀘10

20、.5圖二、四連桿之繪圖其繪出之四連桿為如圖二。黑色為第一桿，藍色為第二桿，紅色為第三桿，綠色 為第四桿。例二、第三桿為結(jié)合桿(coupler)drawlinks(3 2 4 2,0,60,-1,1)圖三、以結(jié)合桿為驅(qū)動桿(r =3 2 4 2)圖三即為所得之答案，此時四連桿為分支型(branch),因為目前之情況無法轉(zhuǎn) 為閉合型，即使將sigma值變號，仍為分支型，如：elf; drawlinks(3 2 4 2,0,60,1,1)圖四、當(dāng)sigma=1時並以結(jié)合桿為驅(qū)動桿(r = 3 2 4 2)利用drawlinks亦可繪出各種角度之圖型，可以作為四連桿運動過程之觀察, 相當(dāng)方便，例如：c

21、lf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,0);endQCPC圖五、多重位置之四連桿運動情形(r =1 3 3 4)其結(jié)果如圖五，當(dāng)然若以第三桿為驅(qū)動桿時，亦可獲得同樣的結(jié)果，例如:clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,-1,1);end圖六、以第三桿為驅(qū)動桿之情形(r = 1 3 3 4)當(dāng)sigma變號時，亦可看出其不同的轉(zhuǎn)動方式，如下：clf;for i=10:20:360,drawlinks(1 3 3 4,0,i,1,1);end圖七、當(dāng)sigma變號(為正時)(r = 1 3 3 4)dra

22、wlinks程式內(nèi)容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 vector for four links%th1: frame angle%th2: crank angle or couple angle%sigma: assembly mode%driver: 0 for

23、 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=1plot(0 rx(1),0 ry(1),'k-','LineWidth',4);hold on;if driver=0plot(0 rx(2),0 ry(2),'b-','LineWidth',1.5); plot(rx(2) rx(3),ry(2) ry

24、(3),'r-','LineWidth',2);elseplot(0 rx(2),0 ry(2),'r-','LineWidth',2); plot(rx(2) rx(3),ry(2) ry(3),'b-','LineWidth',1.5); endplot(rx(1) rx(3),ry(1) ry(3),'-g');plot(rx,ry,'bo');text(0,0,' O');text(rx(1),ry(1),' R');text(

25、rx(2),ry(2),' P');text(rx(3),ry(3),' Q');elsefprintf('Combination of links fail at degrees %6.1fn',th2); endaxis equalgrid on三、drawlimits 函數(shù)四連桿之迴轉(zhuǎn)過程，能完全迴轉(zhuǎn)的情況仍然很少，有些時候無法獲得完整一圈的 迴轉(zhuǎn)。亦即依葛列夫定理四連桿之第一或第二類類型決定，前者為完整迴轉(zhuǎn)型， 後者則有迴轉(zhuǎn)角度之限制，這些限制因四連桿長度決定之。四連桿迴轉(zhuǎn)過程中，有可能其中兩桿會連成一線，或重疊成一線，前者若成立時, 即

26、變成三角形，後者若重疊時，亦會構(gòu)成另一個三角形。理論上連桿構(gòu)成三角形 應(yīng)不會有相對運動。故可稱為四連桿之運動極限。由這兩個極端位置，可以知道 四連桿之最終運動限制。A. 第二桿為驅(qū)動桿在數(shù)學(xué)上，表示這兩個狀況之方法可以利用下列二種不等式進行測試：r1+r2<r3+r4|r1-r2|>|r3-r4|而由其不等式之方向，可以構(gòu)成四種狀況，並進而求得該狀況之角度。下面為第 二桿為驅(qū)動桿時之四種情況：(1)當(dāng) r1+r2三r3+r4, |r1-r2|三 |r3-r4| 時0 min =0 ， 0 max =2 nU)-12 2 2min = - cos r1 +r2 -(r3+r4) /(

27、2 r1 r2),0 max =-12 2 2cos r1 +r2 -(r3+r4) /(2 r1 r2)|h-r2>|r-ci|(3) 當(dāng) r1+r2三r3+r4, |r1-r2|三 |r3-r4| 時0 min =-12 2 2cos r1 +r2 -(r3-r4) /(2 r1 r2),0 max =-12 2 2cos r1 +r2 -(r3+r4) /(2 r1 r2)n +Q>rs + I4(4) 當(dāng) r1+r2三r3+r4, |r1-r2|三 |r3-r4| 時-12 2 20 min = cos r1 +r2 -(r3-r4) /(2 r1 ,-12 2 20 ma

28、x = 2 n - cos r1 +r2 -(r3-r4) /(2 r1 r2)0tJL.dn + i>< c+r-st |B. 第三桿爲(wèi)驅(qū)動桿第三桿結(jié)合桿為驅(qū)動桿時，則仍然取決於四連桿屬於葛列斯荷（Grashof）型或二 型。若屬一型連桿，則當(dāng)?shù)谌龡Ur3為最短桿時，第三桿可以作360度迴轉(zhuǎn)。其 餘之限制條件雖不如以第二桿為驅(qū)動桿者，但其極；限狀況是當(dāng)?shù)诙U與第四桿相平行時，變成無法繼續(xù)迴轉(zhuǎn)，除非它是處於平行四邊形。將四種情況依下列二 不等式之情況加以分類，在這些分類中，若兩式均為等號時，則應(yīng)歸屬於第五類：r1+r3<r2+r4|r1-r3|>|r2-r4|(5) 當(dāng)

29、r1+r3三r2+r4, |r1-r3|三 |r2-r4| 時0 min =0 ， 0 max =2 nI(6) 當(dāng) r1+r3三r2+r4, |r1-r3|三 |r2-r4| 時-12 2 20 min = - cos r1 +r3 -(r2+r4) /(2 r1 r3)，0 max =-12 2 2cos r1 +r3 -(r2+r4) /(2 r1 r3)-12 2 20 min = cos r1 +r3 -(r2-r4) /(2 r1 r3),0 max =-12 2 2cos r1 +r3 -(r2+r4) /(2 r1 r3)(8)當(dāng) r1+r3三r2+r4, |r1-r3|三 |

30、r2-r4| 時-12 2 20 min = cos r1 +r3 -(r2-r4) /(2 r1 ,-12 2 20 max = 2 n - cos r1 +r3 -(r2-r4) /(2 r1 r3)C. fb_a ngle_limits 函數(shù)觀察上面討論之四個極限角度，可以寫一組程式進行計算。由於以第三桿驅(qū)動與 第二桿驅(qū)動，在計算上僅是將其中之r2與r3之位置對調(diào)即可。為尋找上述極限角度0 min、0 max，可用函數(shù)fb_angle_limits進行尋找，其格式如下：function Qstart, Qstop=fb_angle_limits(r,Q1,driver)其中輸入項目有：=

31、四連桿之長度向量，其定義與前函數(shù)相同。Q1 =第一桿之夾角，角度表示（deg）。 driver = 驅(qū)動模式（=0第二桿驅(qū)動；=1第三桿驅(qū)動）。而輸出項為兩個角度：Qstart =驅(qū)動桿（第二桿或第三桿）之最低限角度（deg） Qstop =驅(qū)動桿（第二桿或第三桿）之最高限角度（deg）drawlimits函數(shù)則是呼叫fb_angle_limits.m函數(shù)，然後將其極限位置繪出。其輸 入項目與drawlinks函數(shù)相同。D. 執(zhí)行例（第二桿為驅(qū)動桿） drawlimits（1 2 3 4,0,-1,0）Qstart =3.6e-005Qstop =360圖八、第二桿驅(qū)動，屬全迴轉(zhuǎn)之情形（r =

32、 1 2 3 4）本例中，無論第二桿或第四桿，均可完全迴轉(zhuǎn)。 drawlimits(3 5 2 1,0,1,0)Qstart =-33.557Qstop =33.557圖九、第二類型，右邊有限制，左邊無限(r =3 5 2 1)圖九則因3+5三2+1, |3-5|三|2-1|,故屬前述之第二型限制角度。 drawlimits(5 4 1 3,0,1,0)Qstart =22.332Qstop =51.3183.50.5A/r忌2=51.3p/k/2.51.51.5R0.52.53.54.5圖十、第三類型，前後角度限制(r = 5 4 1 3)本例中，因為各桿長度符合第三類型，即5+4三1+3,

33、 |5-4戶|1-3|,故其角度限制如圖十。drawlimits(4 3 3 5,0,1,0)Qstart =28.955Qstop =331.04圖一、第四類型(r2驅(qū)動)(r =4 3 3 5)此類型與第一類型類似，但角度限制在右邊，因為各桿之長度符合該項條件：4+3 三 3+5, |4-3戸 |3-5。如圖 一。E. 執(zhí)行例(第三桿為驅(qū)動桿)(5) drawlimits(4 5 3 5,0,-1,1)Qstart =3.6e-005Qstop =360圖十二、第三桿與第一桿平行，但可完全迴轉(zhuǎn) (r =4 5 3 5)(6)drawlimits(4 3 3 3,0,-1,1)Qstart

34、=Qstop =-117.28117.28圖十三、第三桿之限制角度相反，屬第六情況(r = 4 3 3 3 )(7) drawlimits(4 4 6 1,0,-1,1)Qstart =26.384Qstop =55.771:q氓J/ Q1/l=26.4Ilc/t S2=55.8/U70.50-0.5-1-1.5-2-2.5-3-3.5-2圖十四、不同角度但同一象限內(nèi)之限制，屬第七狀況(r =4 4 6 1)(8)drawlimits(3 4 4 6,0,-1,1)Qstart =28.955Qstop =331.04-543210-1 -2-3-4圖十五、上下象限之限制角度，兩者和為360度

35、(r =3 4 4 6 )E.程式內(nèi)容functionQstart, Qstop=fb_angle_limits(r,Q1,driver)%function Qstart, Qstop=fb_angle_limits(r,Q1,driver)% Function file that calculates the initial and final angle for four-bar% driver=0 when the crank is the driver.% driver=1 when coupler is the driver.% Variables%r=linkage row vec

36、tor (cm)%Q1=frame angle(degree);%Qstart=initial crank(or coupler)angle (deg)%Qstop=final crank(or coupler)angle (deg)%Programif driver=1, r=r(1) r(3) r(2) r(4);endr仁 r(1);r2=r(2);r3=r(3);r4=r(4);rmin=min(r);rmax=max(r);rtotal=sum(r);if (r1+r2)>(r3+r4 )& abs(r1-r2)<abs(r3-r4) Qstart=acos(-(

37、r4-r3)A2+r2A2+r1A2)/(2*r1*r2); Qstop=acos(r2A2-(r4+r3)A2+r1A2)/(2*r1*r2);endif (r1+r2)>(r3+r4 )& abs(r1-r2)>=abs(r3-r4) Qstart=acos(r2A2-(r4+r3)A2+r1A2)/(2*r1*r2); Qstop=-Qstart;endif (r1+r2)<=(r3+r4)&abs(r1-r2)>=abs(r3-r4)Qstart=0;Qstop=2*pi;endif (r1+r2)<=(r3+r4)&abs(r1-

38、r2)<abs(r3-r4)Qstart=acos(r2A2-(r4-r3)A2+r1A2)/(2*r1*r2);Qstop=2*pi-Qstart;endrange = Qstop-Qstart;Qstart=Q1+(Qstart+0.0000001*range)*180/pi;Qstop=Q1+(Qstop-0.0000001*range)*180/pi;if Qstop<Qstart;temp=Qstart;Qstart=Qstop;Qstop=temp;end% %function drawlimits(r,th1,sigma,driver)%function drawlm

39、its(r,th1,sigma,driver)%draw the positions of four-bar links%call f4bar funcion%r:row vector for four links%th1: frame angle%sigma: assembly mode %driver: 0 for crank, 1 for couplerQstart, Qstop=fb_angle_limits(r,th1,driver)clf;values b=f4bar(r,th1,Qstart,0,0,sigma,driver);rr=values(:,1); rr(3)=rr(1)+rr(4);rx=real(rr);rx(4)=0; ry=imag(rr);ry(4)=0;if b=1plot(0 rx(1),0 ry(1),'k-','LineW